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The  paper  reports  on  a  theoretical  investigation  of  the  conver¬ 
gence  properties  of  several  finite  element  approximations  in  current 
use  and  assesses  the  magnitude  of  the  principal  errors  resulting  from 
their  use  for  certain  classes  of  structural  problems.  The  method  is 
based  on  classical  order  of  error  analyses  commonly  used  to  evaluate 
finite  difference  methods.  Through  the  use  of  the  Taylor  series  dif¬ 
ferential  or  partial  differential  equations  are  found  which  represent 
the  convergence  and  principal  error  characteristics  of  the  finite 
element  equations.  These  resulting  equations  are  then  compared  with 
known  equations  governing  the  continuum,  and  the  error  terms  are 
evaluated  for  selected  problems.  Finite  elements  for  bar, beam,  plane 
stress,  and  plate  bending  problems  are  studied  as  well  as  the  use  of 
straight  or  curved  elements  to  approximate  curved  beams.  The  results 
of  the  study  provide  basic  information  on  the  effect  of  interelement 
compatibility,  unequal  size  elements,  discrepancies  in  triangular 
element  approximations,  flat  element  approximations  to  curved 
structures,  and  the  number  of  elements  required  for  a  desired  degree 
of  accuracy. 


995 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

OCT  1968 

2.  REPORT  TYPE 

3.  DATES  COVERED 

00-00-1968  to  00-00-1968 

4.  TITLE  AND  SUBTITLE 

5a.  CONTRACT  NUMBER 

Accuracy  and  Convergence  of  Finite  Element  Approximations 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Flight  Dynamics  Laboratory, Wright  Patterson  AFB, OH, 45433 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

See  also  AD0703685,  Proceedings  of  the  Conference  on  Matrix  Methods  in  Structural  Mechanics  (2nd) 

Held  at  Wright-Patterson  Air  Force  Base,  Ohio,  on  15-17  October  1968. 

14.  ABSTRACT 

15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 
ABSTRACT 

18.  NUMBER 
OF  PAGES 

34 

19a.  NAME  OF 
RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


AFFDL-TR-68-150 


SECTION  I 
INTRODUCTION 


Finite  element  methods  have  been  used  for  many  years  with  success  in  the  analysis  of 
complex  structures  and  many  aerospace  structures  are  designed  on  the  basis  of  these  analyses. 
In  spite  of  such  widespread  use,  not  enough  is  known  of  the  theoretical  accuracy  and  conver¬ 
gence  properties  of  these  models  when  used  to  represent  a  structure.  Accuracy  studies  are 
usually  based  on  numerical  solutions  to  restricted  problems  for  comparison  with  known 
results.  Convergence  studies  are  carried  out  by  investigating  the  convergence  of  the  numerical 
results  as  the  number  of  elements  is  increased  (Reference  1).  Such  methods,  while  valuable  in 
providing  a  cursory  assessment  of  the  reliability  of  various  model  approximations,  are  heavily 
dependent  on  the  numerical  data  and  the  problem  studied  and  may  obscure  the  true  character 
of  the  approximation.  More  basic  information  is  needed  on  the  theoretical  convergence 
properties  of  various  finite  elements  for  structural  approximations. 

The  purpose  of  the  present  paper  is  to  report  on  a  theoretical  investigation  of  the  accuracy 
of  several  stiffness  finite  element  method  approximations  in  current  use  and  to  assess  the 
magnitude  of  the  errors  resulting  from  the  use  of  these  approximations  for  certain  classes  of 
structural  problems.  The  present  paper  documents  preliminary  results  on  convergence  of 
several  models  given  orally  in  Reference  2,  extends  the  study  to  additional  models,  and 
utilizes  this  data  for  accuracy  investigations.  The  method  used  is  based  on  classical  order  of 
error  analyses  commonly  used  to  evaluate  the  discretization  errors  of  finite  difference 
methods.  Through  the  use  of  the  Taylor  series,  the  ordinary  or  partial  differential  equations 
are  found  from  which  the  convergence  and  error  characteristics  of  the  finite  element  equations 
can  be  determined.  These  resulting  equations  are  then  compared  with  the  known  equations 
governing  the  structure.  An  estimate  of  the  discretization  error  in  the  finite  element  approxi¬ 
mation  is  evaluated  by  procedures  given  in  References  3  and  4  for  a  limited  class  of  deflection 
and  vibration  problems  to  provide  simple  formulas  for  the  size  of  an  element  required  to 
obtain  a  certain  degree  of  accuracy.  Finite  elements  for  bar,  beam,  plane  stress,  and  plate 
bending  problems  are  studied  as  well  as  the  use  of  straight  and  curved  elements  to  approxi¬ 
mate  a  curved  beam. 
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SECTION  n 

ERROR  ANALYSIS  PROCEDURE 

Two  major  sources  of  error  result  from  the  use  of  finite  element  methods  to  solve 
structural  problems.  These  may  be  conveniently  separated  into  round-off  error  and  dis¬ 
cretization  error.  Round-off  error  is  that  error  associated  with  the  accuracy  with  which 
numbers  are  manipulated  in  a  computer  and  is  not  considered  in  the  present  paper.  Dis¬ 
cretization  error  is  that  error  associated  with  using  discrete  variables  to  represent  a 
problem  where  the  state  variables  are  continuous.  This  error  occurs  irrespective  of  the 
accuracy  of  numerical  calculations  and  occurs  in  structural  problems  when  finite  elements 
are  used  to  approximate  a  continuous  structure.  Discretization  error  may  be  of  two  kinds: 
(1)  errors  proportional  to  the  size  of  the  element  which  vanish  as  the  element  size  vanishes 
and  (2)  errors  which  do  not  vanish  when  the  element  size  vanishes.  Elements  and/or  patterns 
which  lead  to  the  second  kind  of  discretization  errors  are  unsatisfactory  approximations 
and  should  be  recognized  and  avoided.  The  present  paper  deals  with  an  assessment  of 
discretization  error  in  finite  element  approximations. 

The  method  used  in  the  study  is  to  obtain  the  typical  finite  element  equations  which 
express  force  equilibrium  at  a  reference  node  point  in  terms  of  displacement  variables. 
These  finite  element  equations  (which  are  a  class  of  difference  equations)  are  then  expanded 
in  Taylor  series  about  the  nodal  point  to  obtain  the  differential  equations  equivalent  to  the 
finite  element  equations  at  that  node.  The  resulting  differential  equations  are  compared  with 
the  governing  equations  for  the  continuum  approximated.  A  simple  bar  element  approximation 
is  treated  in  detail  as  an  example  to  characterize  the  method  and  to  define  the  terms  to  be 
used  in  the  study. 

DISCRETIZATION  ERRORS 


The  force-displacement  relations  of  a  typical  one -dimensional  structural  element  having 


ends  i-1  and  i  are 


‘  Fi 

_  Fi 


8; 


where  F  .  and  8 .  are  the  vector  of  nodal  forces  and  displacements  at  the  ith  node  and  K  is 
the  element  stiffness  matrix.  Consider  a  bar  of  constant  cross  sectional  area  A  subjected 
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to  a  distributed  axial  load  p(x)  and  approximated  by  finite  elements,  where  x  denotes  distance 
along  the  bar  (Figure  1).  For  a  bar  finite  element  F.  is  the  nodal  extensional  force,  8.  is  the 
corresponding  displacement  u.  and 

«  =¥-[-!  'I  ] 

where  h  is  the  length  of  the  element  and  E  is  Young’s  modulus. 


In  finite  element  methods  distributed  loads  are  replaced  by  concentrated  loads  at  the 
node  points.  There  are  several  rational  ways  in  which  distributed  loads  may  be  converted  to 
concentrated  loads.  The  primary  purpose  of  the  present  paper  is  to  investigate  the  accuracy 
of  finite  element  approximations  to  structures  and  this  can  be  done  by  using  a  very  simple 
lumping  procedure  for  loads.  This  procedure,  for  example,  for  the  bar  gives  the  concentrated 
load  as  the  value  of  the  distributed  load  at  a  grid  point  multiplied  by  one-half  of  the  total 
length  of  the  two  adjoining  elements.  Treating  the  distributed  load  in  this  manner  and  utilizing 
Equations  1  and  2  leads  to  the  following  equation  for  equilibrium  at  a  typical  i  point  located 
between  two  segments  of  length  h  and  ah  (Figure  1). 
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Equation  3  is  a  typical  finite  element  equilibrium  equation  for  the  indicated  approximation. 


The  behavior  of  the  system  of  Equation  3  is  investigated  as  the  number  of  equations  goes 
to  infinity  and  the  size  of  the  element  vanishes.  This  is  done  by  examining  Equation  3  in  the 
limit  as  h  approaches  zero  with  the  aid  of  the  Taylor  series  expansion  of  displacements  at 
points  i-1  and  i+1  about  the  ith  point.  This  procedure  results  in 
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where  primes  denote  differentiation  with  respect  to  x  and  where  the  subscript  i  has  been 
omitted  from  u..  (Omission  of  such  subscripts  is  done  consistently  throughout  the  paper). 
Equation  4  is  the  differential  equation  equivalent  to  the  finite  element  equation  at  the  i  node. 
The  terms  in  Equation  4  which  are  not  multiplied  by  powers  of  h  comprise  exactly  the 
governing  differential  equation  for  the  continuous  bar  and  the  remaining  terms  are  the  dis¬ 
cretization  errors  resulting  from  use  of  the  finite  element  approximations. 


If  a  finite  element  equation  converges  to  the  governing  differential  equation  for  the  con¬ 
tinuous  structure  as  h  vanishes,  the  finite  element  approximation  will  be  defined  in  this  paper 
to  be  a  consistent  approximation.  Those  terms  in  Equation  4  which  differ  from  the  continuous 
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P(x) 

(a)  Continuous  "bar. 
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ah 


i+1 
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(b)  Approximation  of  continuous  bar  by  finite  elements. 


F1'V 


(c)  Forces  and  displacements  m  bar  finite  elements. 


[ 


(d) 


Equilibrium  at  i^ 


grid  point  in  finite  element  assemblage. 


Figure  1.  One  Dimensional  Bar  and  Its  Finite  Element  Representation 

structure  equation  are  the  discretization  errors  and  the  principal  error  in  the  approximation 
is  the  set  of  terms  in  the  discretization  error  multiplied  by  the  lowest  power  of  h.  The  power 
of  h  multiplying  the  principal  error  denotes  the  order  of  the  discretization  error  of  the  finite 
element  equation  and  the  rate  of  convergence  as  h  vanishes.  Thus,  the  finite  element  equation 
leading  to  equation  (4)  has  a  discretization  error  of  order  h.  Note  that  if  the  segments  are 
equal  (  a  =  1),  the  discretization  error  is  of  order  h2.  Similarly,  if  the  finite  element  equations 
do  not  converge  to  the  governing  differential  equations,  the  discretization  error  would  be  of 
order  h^  or  one. 

HARMONIC  LOADING  AND  VIBRATION  EXAMPLE 

The  error  analysis  procedure  described  in  the  previous  section  gives  only  the  rate  of 
convergence  of  the  finite  element  approximation.  The  magnitude  of  the  discretization  error 
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in  the  bar  finite  element  approximation  can  be  evaluated  for  the  special  case  of  a  bar  supported 
at  each  end,  subjected  to  a  sinusoidally  distributed  static  loading 

m  7T  x 


p  =  p  si  n  - 

o  L 


(5) 


and  approximated  by  equal  length  elements.  In  Equation  5  pQ  is  the  amplitude  of  the  loading, 
m  the  number  of  half  waves  and  L  the  length  of  the  bar.  For  this  case  Equation  4  becomes 
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If  the  loading  is  regarded  as  a  continuous  function  as  was  done  in  Reference  3  and  4,  the 
solution  to  Equation  6  is 

m-rrx 

u  :  u  sin - 

o  L 

where  uq  is  the  amplitude  of  the  displacement.  Substitution  of  Equation  7  into  Equation  6 
results  in 
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where 
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is  the  exact  solution  to  the  bar  for  the  sine  loading, 
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is  the  principal  error  in  deflection  resulting  due  to  the  finite  element  approximation  and 

N  -  L/m  (ID 

N  h 

is  the  number  of  elements  per  harmonic  half-wave  used  to  approximate  the  bar.  The  following 
are  results  for  various  values  of  N  and  error 


N 

3 

4 
9 


€ 

0.10 

.05 

.01 


1000 


AFFDL-TR-68-150 


Thus,  approximately  three  elements  per  deflection  half-wave  are  needed  to  keep  the  error 
in  finite  element  deflection  calculations  at  a  node  to  within  ten  percent. 


A  similar  approach  can  be  used  to  determine  the  error  in  natural  frequency  of  the  bar 
example  when  approximated  by  equal  length  elements.  For  vibration  behavior  where  a  simple 
lumping  process  is  used  to  obtain  a  diagonal  mass  matrix  the  counterpart  of  Equation  6  is 


12 


1  v  * 
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EA 


u  =  0 


(12) 


where  m  is  the  mass  per  unit  length  and  cu  is  the  circular  frequency.  The  vibration  mode 
shape  is  the  same  as  Equation  7  and  the  eigenvalues  of  Equation  12  are 


where 


(13) 


UJ 


ex 
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and  £  is  the  discretization  error  due  to  the  finite  element  approximation  given  by  Equation  10. 
Note  that  the  error  in  the  square  of  frequency  is  the  same  as  that  resulting  from  static  de¬ 
flection  calculations  except  of  opposite  sign.  Thus,  the  frequency  calculations  converge  from 
below  and  deflection  calculations  converge  from  above. 


For  completeness,  Table  1(a)  summarizes  the  principal  error  terms  for  the  general  bar 
finite  element  approximations  and  Table  II  evaluates  these  errors  for  the  harmonic  response 
examples. 
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SECTION  HI 

RESULTS  AND  DISCUSSION 


The  method  described  in  the  previous  section  was  used  to  investigate  the  convergence  of 
beam  elements,  straight  and  curved  element  approximations  to  an  arch,  rectangular  and 
triangular  plane  stress  elements,  and  plate  bending  elements.  The  results  of  the  investigation 
are  summarized  in  this  section  together  with  a  general  discussion.  A  summary  of  the  error 
study  details  is  given  in  Tables  I  and  II. 

ONE  DIMENSIONAL  BENDING  ELEMENT 


When  bending  behavior  is  introduced  finite  element  models  and  the  discretization  error 
analysis  procedure  become  more  complex.  To  indicate  the  additional  features  brought  on  by 
bending,  consider  a  simple  prismatic  beam  subjected  to  lateral  pressure  q  and  approximated 
by  an  assemblage  of  beam  bending  finite  elements  of  length  h  (Figure  2).  The  finite  element 
nodal  variables  for  this  problem  are 


».  ■  [  e,  } 

where  w.  and  8 .  are  the  displacement  and  rotation  at  a  node,  and  two  finite  element  equilib¬ 
rium  equations  are  obtained  at  each  node.  These  two  equations  are  expanded  in  a  Taylor  series 
about  the  ith  point  (considering  both  w  and  6  as  independent  variables)  giving 
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where  I  is  the  moment  of  inertia  of  the  beam  cross  section.  Since  beam  behavior  for  a  con¬ 
tinuum  is  defined  by  only  one  independent  variable  w,  it  is  useful  to  eliminate  8  in  so  far  as 
possible  from  the  finite  element  equation.  This  is  done  by  differentiating  Equation  17  to  obtain 
expressions  for  derivatives  of  8  and  sequentially  back  substituting  these  derivatives  into 
both  Equations  16  and  17. 
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TABLE  I.-  FINITE  ELEMENT  DISCRETIZATION  ERROR 


(a)  One-Dimensional  Elements 


Element 


Nodal 

varieties 


Governing 

differential  equations 


Error  terms 

(appearing  in  left-hand 
side  of  equations) 


(a)  Straight  segments 
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TABLE  I.-  UNITE  ELEMENT  DISCRETIZATION  ERROR  -  Continued 
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TABLE  I.-  FINITE  ELEMENT  DISCRETIZATION  ERROR  -  Concluded 
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TABLE  II.-  ERROR  TERMS  FOR  HARMONIC  RESPONSE  EXAMPLES 

Element  Error  terms 
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Figure  2.  Beam  Element 
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Equations  16  and  17  finally  can  be  put  in  the  form 
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(19) 


Equations  18  and  19  show  that  in  the  limit  as  h  vanishes  the  beam  finite  element  equations 
converge  to  the  familiar  beam  equation  and  the  correct  constraint  equation  between  the  two 
finite  element  variables  6  and  w.  The  equations  also  show  that  the  principal  error  term  in 
the  approximation  is  order  h4  (Table  1(a)). 


The  principal  error  was  evaluated  for  a  finite  element  approximation  to  a  simply  sup¬ 
ported  beam  of  length  L  subjected  to  a  sinusoidally  distributed  lateral  load  q. 

m  7TX 

q  1  q0  sin  — —  (20) 

where  qQ  is  a  constant.  The  error  in  the  lateral  deflection  w  of  the  beam  due  to  the  finite 
element  approximation  is  (Table  II). 
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Equation  21  also  gives  the  magnitude  of  the  error  in  frequency  determination  due  to  the  finite 
element  approximation  (Table  II). 

APPROXIMATION  OF  CURVED  STRUCTURES 


Straight  elements  are  often  being  used  to  approximate  curved  structures  such  as  curved 

beams,  arches,  and  shells.  To  gain  some  insight  into  the  influence  of  curvature,  an  arch  of 

radius  R  was  approximated  by  conventional  straight  beam  elements  which  also  had  extensional 

capability.  The  arch  loading  is  a  normal  pressure  q  and  a  tangential  distributed  loading  p 

til 

(Figure  3).  At  a  typical  i  node  three  variables  are  required  to  define  element  behavior. 
For  this  problem  it  is  convenient  to  take  these  variables  as 
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Figure  3.  Arch  Approximately  by  Straight  and  Curved  Elements 
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where  u  and  w  are  the  tangential  and  radial  displacements  and  8  the  rotation.  Three  finite 
element  equations  result  from  force  and  moment  equilibrium.  As  the  element  size  vanishes, 
the  moment  equilibrium  equation  at  the  i^1  node  converges  to  the  correct  constraint  equation 
between  the  rotation  8,  radial  displacement  w,  and  tangential  displacement  u 


Mw'+-s-)=o 


(23) 


and  the  tangential  and  normal  equilibrium  equations  converge,  respectively,  to 
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In  Equations  24  and  25  the  rotation  8  has  been  eliminated  by  using  Equation  23  in  a  manner 
similar  to  that  done  earlier  for  the  other  bending  problems.  A  typical  set  of  arch  continuum 
equations  such  as  those  given  in  Reference  5  is  composed  of  only  the  left-hand  side  of 
Equations  23,  24,  and  25;  namely,  equation  23  plus 


E  A  (-  u"  +■  ~  )  -  p  s  0 


(26) 


El(‘iv  +  2^  +  ^)+“(-f +  ^)--0 

Comparison  of  Equations  24  and  25  with  26  and  27  shows  that  finite  element  equations  re¬ 
sulting  from  a  straight  element  approximation  do  not  converge  to  these  arch  equations  since 
the  right  side  of  24  and  25  does  not  vanish.  Thus,  the  discretization  error  appears  to  be  of 
order  one  and  suggests  that  some  errors  result  from  the  use  of  straight  elements  to  approx¬ 
imate  the  bending  behavior  of  the  curved  structure. 


A  study  was  also  made  of  the  use  of  curved  elements  to  approximate  the  curved  structure. 
While  the  details  have  been  omitted  here  the  stiffness  matrix  for  a  curved  element  was 
derived  from  the  strain  energy  consistent  with  arch  Equations  26  and  27  (Reference  5).  The 
displacements  were  approximated  by  assuming  that  arch  tangential  and  normal  displacements, 
were  linear  and  cubic,  respectively,  over  the  curved  element  length.  The  resulting  finite 
element  equations  were  investigated  and  it  was  found  that  the  element  pattern  converged  to 
the  arch  Equations  26  and  27  with  an  error  of  order  h2. 


1011 


AFFDL-TR-68-150 


An  assessment  of  the  magnitude  of  the  order  one  discretization  error  terms  in  Equations 
24  and  25  can  be  obtained  by  considering  a  closed  circular  ring  subjected  to  a  harmonic 
lateral  loading  and  approximated  by  straight  elements.  In  the  limit  as  the  element  size  h 
vanishes  the  resulting  straight  element  psuedo  arch  is  governed  by  Equations  24  and  25  and 

p  =  o  (28) 

q  =  qo  Sin  {29) 

The  principal  errors  in  the  u  and  w  deflections  and  ew  due  to  the  finite  element  approxi¬ 
mation  are  (Table  II) 


where  p  is  the  radius  of  gyration  of  the  arch.  Equations  28  through  31  indicate  that  the  order 

one  discretization  errors  lead  to  errors  in  deflection  which  are  proportional  to  the  square  of 

the  ratio  of  the  radius  of  gyration  to  the  arch  radius.  For  a  typical  arch  cross  section  these 

errors  are  quite  small  and  well  within  accuracy  requirements  for  engineering  purposes 

provided  enough  elements  are  used.  Tangential  displacement  errors  are  also  proportional 
o 

to  m  ;  however,  the  ring  theory  breaks  down  for  very  high  harmonics. 

Since  the  magnitude  of  the  order  one  errors  are  proportional  to  ^  j  and  therefore 
quite  small  for  thin  arches  it  suggests  that  Equations  24  and  25  are  acceptable  arch  equations. 
Further  investigation  shows  that  a  suitable  thin  arch  theory  can,  in  fact,  be  derived  which  leads 
to  Equations  24  and  25.  This  theory  canbe  obtained  from  the  arch  theory  of  Reference  5  if  the 
change  in  curvature  are  modified  by  an  additional  term  which  is  the  extensional  strain 
divided  by  the  arch  radius.*  According  to  the  Koiter  criteria  for  thin  shells  (Reference  12) 
such  modifications  are  admissible  variations  of  a  first  approximation  theory  of  thin  shells 
and  it  seems  reasonable  to  apply  this  criteria  to  arch  theory.  Thus  both  straight  and  curved 
elements  provide  a  convergent  approximation  to  a  first  approximation  thin  arch  theory  as 
the  element  size  vanishes. 


*The  authors  are  indebted  to  Professor  B.  Budiansky  of  Harvard  University  for  suggesting 
that  alternate  first  approximation  arch  theories  be  investigated. 
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TWO  DIMENSIONAL  PLANE  STRESS  ELEMENTS 

The  linear  elastic  plane  stress  equations  for  equilibrium  in  the  x  and  y  directions 
formulated  in  terms  of  displacements  are,  respectively: 


u’xx  u. 

,  1  +  M  .  px 

'yy*  2  V’*y+B  "° 

(32) 

1  4-  /i 

2  U’xy 

'"M  ,  ,  py 

2  'X1<  fv’yyf  B  ‘° 

(33) 

Here  u  and  v  are  the  displacements  in  the  x  and  y  directions,  respectively,  p  ,  p  the  dis- 

/  o\  x  y 

tributed  forces,  H-  Poisson’s  ratio,  B  =  Et/(l  -fx  )  the  extensional  stiffness  and  t  the  plate 
thickness.  A  satisfactory  finite  element  approximation  should  lead  to  equations  which  converge 
to  Equations  32  and  33  at  a  node  as  the  element  size  vanishes. 

Rectangular  Elements 

Two  rectangular  plane  stress  plate  elements  were  investigated  in  the  study  and  their 
stiffness  properties  are  documented  in  reference  6.  For  the  first  model,  denoted  linear  stress 
model,  the  stresses  in  the  x  and  y  directions  are  assumed  to  vary  linearly  while  the  shear 
stress  is  constant  (Figure  4).  For  the  second  model,  denoted  a  linear  edge  displacement 
model,  the  displacements  along  an  edge  of  the  element  are  assumed  to  vary  linearly.  The 
nodal  variables  used  to  define  the  stiffness  matrices  for  these  finite  elements  are 


and  a  typical  finite  element  equation  contains  contributions  from  all  elements  contiguous  to  the 

node.  The  pattern  arrangement  composed  of  equal  elements  is  shown  in  Figure  4.  The 

distributed  forces  on  the  plane  stress  body  were  again  concentrated  in  a  simple  fashion  based 

on  the  value  of  the  distributed  force  at  each  node  location.  The  typical  finite  element  equations 

were  obtained  and  the  error  terms  evaluated.  An  investigation  of  the  convergence  of  the 

finite  element  equilibrium  equation  for  both  models  showed  them  to  converge  to  the  plane 

o 

stress  equations  32  and  33  with  a  principal  error  of  order  h  .  (Table  1(b)). 


The  principal  error  was  evaluated  for  the  case  where  the  two  rectangular  elements 
were  used  to  approximate  a  square  plate  in  plane  stress  subjected  to  a  harmonic  in-plane 
loading,  in  the  y  direction.  The  loading  is 

px  z  0  (34) 


rmrx  mr 
py  =  po  sm  a~  sir>  ~ 
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and  the  plate  is  supported  on  the  boundary  such  that  the  force  resultant  in  the  x  direction  and 

the  v  displacement  both  vanish.  For  /a  =  0.3  and  m  =  n  the  errors  in  the  u  and  v  displace¬ 

ments,  ffu  and  respectively,  are  (Table  II). 

1.  linear  stress  model 

1.85 

€u  =  ~ (36) 
2.40 

'  N2  (37) 

2.  linear  edge  displacement  model 

I.  15 

€u  (38) 

I.  97 

€v=  n2  (39) 

where  N  is  the  number  of  elements  per  half  wave-length. 

Triangular  Elements 


Results  were  also  obtained  for  the  use  of  the  classical  triangular  plate  element  (Ref¬ 
erence  7)  to  approximate  plane  stress  problems.  Arrangements  or  patterns  A,  B,  anri  c 
were  investigated  for  the  convergence  of  three  right  triangular  elements  (Figure  5  for 


patterns  about  ij  point).  It  was  found  that  pattern  A  converges  to  the  required  plane  stress 
equations  32  and  33  and  the  principal  error  is  of  order  h2.  On  the  other  hand,  the  corre¬ 
sponding  equations  for  pattern  B  are 


u 


xx 


I  -  H- 
2 


’yy 


(40) 


1  -  /A  Py 

2  V,xx  +  v>yy  +  =  0  (4() 


and  pattern  C  are 


xx 


I  ~  A6 
2 


yy 


+  ( i  V, 


xy 


=0 


(42) 


(l  +/A)  u ’xy  +  2  ^  v,xx  +  v’yy  +  =  0 


The  additional  error  terms  for  all  patterns  proportional  to  h2  are  given  in  Table 


(43) 

1(b). 
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Some  points  can  be  noted  by  comparing  the  convergence  characteristics  when  h  =  0  of 
the  B  and  C  equations  with  Equations  32  and  33.  First  of  all  both  patterns  B  and  C  lead  to  a 
discretization  error  of  order  one.  The  pattern  B  equations  do  not  contain  cross -derivative 
terms  suggesting  that  convergence  of  shear  behavior  at  the  nodal  point  is  poor.  This  is  not 
unexpected  since  there  is  no  mechanism  in  the  finite  element  equations  for  representing 
changes  in  shear  at  the  nodal  point  due  to  the  arrangement  of  the  elements.  On  the  other 
hand,  the  pattern  C  equations  over  prescribe  the  cross-derivative  term  by  a  factor  of  two. 
Note  that  the  difficulty  arises  from  the  element  arrangement  rather  than  the  element 
properties  since  the  element  used  here  fully  represents  all  states  of  plane  stress.  This 
indicates  that  convergence  difficulties  may  arise  for  some  triangular  elements  as  a  result  of 
poor  element  arrangement  even  though  the  element  is  well  formulated. 

Since  the  difficulty  with  pattern  B  is  due  to  its  inability  to  represent  the  cross  derivative 
at  the  node,  better  convergence  properties  would  be  expected  if  additional  degrees  of  freedom 
were  used  to  characterize  the  right  triangular  element  behavior.  Such  added  degrees  of 
freedom  might  be  the  deflections  at  the  midpoint  of  the  various  edges  or  the  derivatives  of 
displacements  at  nodes. 

Fortunately,  if  patterns  B  and  C  are  used  in  structural  idealizations,  they  usually  occur 
in  pairs  (Figure  5)  and  the  under  prediction  of  the  shear  stiffness  at  one  point  is  compensated 
for  to  some  extent  by  an  over  prediction  of  shear  stiffness  at  a  neighboring  point.  Nevertheless, 
these  results  do  suggest  that  caution  should  be  exercised  to  ensure  that  an  excessive  number 
of  either  patterns  A  or  B  does  not  occur  in  a  structure  when  the  results  are  strongly  de¬ 
pendent  on  shear  stiffness,  A  more  consistent  approach  is  to  use  pattern  A  since  it  converges 
to  the  appropriate  plane  stress  equation. 

For  completeness,  the  convergence  of  a  pattern  composed  of  equilateral  triangles  in 
plane  stress  was  also  investigated.  The  typical  pattern  is  shown  in  Figure  6  and  the  resulting 
finite  element  equations  were  found  to  converge  to  the  plane  stress  equations  with  a  principal 
error  of  order  h2.  (Table  1(b)). 

TWO-DIMENSIONAL  PLATE  BENDING  ELEMENTS 

Three  rectangular  plate  bending  models  were  investigated  to  determine  the  convergence 
and  principal  errors  of  the  resulting  finite  element  equations.  The  models  investigated  were 
those  developed  by  Papenfuss  (Reference  8),  Melosh  (Reference  9),  and  one  developed 
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Figure  6.  Equilateral  Triangular  Finite  Element  Pattern 


independently  by  Adini  and  Clough  {Reference  1  or  Reference  10)  and  Melosh  (Reference  11). 
(These  will  be  denoted  respectively  as  the  Papenfuss,  Melosh,  and  ACM  models.)  The  pattern 
arrangement  is  shown  in  Figure  7.  The  nodal  variables  for  these  finite  element  models  are 


where  w  is  the  lateral  displacement  and  9  and  the  rotations  about  the  y  and  x  axes, 
respectively.  On  the  basis  of  the  beam  results  a  consistent  set  of  three  plate  bending  finite 
element  equations  should  be  expected  to  converge  to 

V4  w  =  -SL  (45) 

=  w,  x  (46) 


9 
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as  the  element  size  vanishes.  Here  Equation 45  is  the  familiar  plate  equation  and  Equations  46 
and  47  are  constraint  equations  between  the  rotations  8  and  <fi  and  derivative  of  w.  The  finite 
element  equilibrium  equations  for  the  three  models  at  a  node  point  were  expanded  in  a  Taylor 
series  and  manipulated  to  a  convenient  form  in  a  manner  similar  to  that  done  for  the  beam 
equations.  All  three  models  lead  to  equations  which  converge  to  constaint  equations  of  the 


form  of  Equations  46  and  47  and  the  Melosh  and  ACM  models  also  converge  to  the  Plate 

Equation  45.  The  Papenfuss  model,  however,  converges  to  a  psuedo  plate  equation 

2 


V4w  + 1 


2  a 


35  35a' 


25 )  "’  xxyy 


w 


(48) 


where  a  is  the  aspect  ratio  of  the  element  and  thus  has  a  principal  error  of  order  one.  The 

2 

principal  errors  for  all  other  equations  for  the  three  models  are  proportional  to  h  with  the 
exception  of  the  constraint  equations  for  the  ACM  model  which  are  proportional  to  h4 
(Table  1(c)). 


It  is  well  known  from  numerical  calculations  (Reference  1)  that  the  Papenfuss  model  has 
some  deficiencies  and  that  the  source  of  the  discrepancies  is  the  inability  of  the  model  to 
describe  the  twist  behavior  of  a  plate.  This  discrepancy  term  shows  up  as  an  incorrect  cross¬ 
derivative  term  in  Equation  48  when  h  vanishes. 


Square  elements  of  the  three  models  were  then  used  to  approximate  a  simply  supported 
square  plate  subjected  to  a  harmonic  loading. 


m  7TX  n  7TX 

q  1  «<> sin  sin  T~ 


(49) 


The  same  procedure  was  used  to  approximate  the  lateral  vibration  characteristics  of  the 
plate.  The  error  in  deflection  €  ^  and  the  error  in  frequency  ^ue  the  element 
approximation  is,  for  m  =  n  and  /x  =  0.3  (Table  II). 


Popenf u  ss 

.2646 


,0  463  - 


N* 


-.0486  - 


.2909 

N2 


Melosh 

.7  08 
N  2 

.708 

N2 


ACM 

t.069 

N2 

I  .069 
N2 
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where  N  =  is  the  number  of  elements  per  Fourier  half-wave.  Figure  8  shows  a  sketch 
of  the  plate  and  a  plot  of  the  ratio  of  the  finite  element  results  to  the  exact  result  for  the 
three  elements  as  the  l/N  vanishes.  The  results  show  that  in  the  limit  as  the  size  of  the 
element  vanishes,  the  results  for  this  problem  based  on  the  Papenfuss  model  have  an  error 
of  approximately  five  percent.  Note  also  that  the  Papenfuss  model  converges  from  below 
and  the  other  two  models  converge  from  above. 


Since  the  Papenfuss  model  equations  do  not  converge  to  the  plate  equation,  a  solution  was 
obtained  for  the  resulting  Papenfuss  psuedo  plate  Equation  48  for  a  rectangular  planform 
subjected  to  a  uniform  load  q,  having  simple  support  boundary  conditions  and  approximated 
by  elements  having  the  same  aspect  ratio  as  the  plate.  The  solution  was  obtained  by  classical 
Fourier  series  expansion  for  the  load  q  and  deflection  shape  as 


T  =  I  I 

m  =1  n  =  I 

*  -  I  I 

m  =  I  n  =  I 

where  a  and  b  are  the  lengths  of  the  plate  in 


m7rx  n7ry 

<lmn  Sm  1 —  Si"  r- 

m  n  o  b 


m7rx  riTry 

w  sin  -  sin  — — 

mn  o  b 


the  x  and  y  direction,  respectively. 


(50) 


(51) 


The  w^n  are  related  to  the  q  by 


1 6  q 


w 


mn 


T  0  mn 


r  4  2  Z 

rm  m  n 

/o  ,  2a*  2  ,  2 

!  +  —  1 

L  4  +  2  2 

a  a  b 

1  L  T  f  .  T 

'  35  35a2  25 

1  J 

The  following  compares  the  center  deflection  w  of  the  Papenfuss  psuedo  plate  with  the  exact 

c 

plate  results 

WC  0  3 

C  X  to3 


—  4 

qa 


a  /b  =  I 
a  /to  =  1/2 


Papenfuss  Plate 
3.87 
9-64 


E  xoct 
4.06 
10.13 


These  results  indicate  that  the  Papenfuss  psuedo  plate  has  an  error  of  approximately  five 
percent  for  both  cases.  Numerical  results  obtained  in  Reference  1  for  the  above  cases  appear 
to  be  converging  toward  these  analytical  results. 
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COMPARISON  OF  ELEMENT  MODELS 

Some  general  results  can  be  deduced  by  comparing  the  convergence  and  principal  error 
results  of  the  various  element  approximations.  One  result  deals  with  the  requirement  of 
interelement  compatability.  Consideration  of  the  elements  for  both  plane  stress  and  bending 
gave  examples  where  this  requirement  was  neither  necessary  nor  sufficient  to  insure 
convergence.  The  rectangular  linear  stress  model  is  not  interelement  compatible  in  dis¬ 
placements  and  the  rectangular  Melosh  and  ACM  bending  models  are  not  compatible  in  slope; 
yet  these  three  lead  to  convergent  equations.  On  the  other  hand,  the  right  triangular  plane 
stress  patterns  A  and  B  are  interelement  compatible  in  displacements  and  the  Papenfuss 
plate  bending  element  is  compatible  in  slope  and  displacements  and  yet  these  elements  and 
arrangements  do  not  lead  to  convergent  equations. 

Some  influence  of  unequal  length  segments  is  seen  from  the  change  in  principal  errors 
for  the  bar  approximations.  For  bar  elements  of  equal  length  the  principal  error  is  propor¬ 
tional  to  h  and  for  bars  of  unequal  length  proportional  to  h.  From  the  asymmetric  character 
of  the  Taylor  series  expansion  about  the  reference  point  similar  reduction  of  the  order  of 
error  would  be  expected  for  other  elements.  This  slower  rate  of  convergence  suggests  that 
results  may  be  less  accurate  when  structures  are  approximated  by  unequal  elements  than 
where  approximated  by  equal  length  segments. 

Comparison  of  the  errors  for  the  plane  stress  with  those  of  the  bar  subjected  to  harmonic 
loading  indicates  that  approximately  nine  bar  elements  and  15  square  plane  stress  elements 
per  half  wave  in  one  direction  are  required  for  one  percent  error.  Approximately  two  beam 
elements  and  nine  square  Melosh  plate  bending  elements  are  required  for  1  percent  error. 
This  result  indicates  that  more  elements  are  required  per  wave  length  in  one  direction  for 
two-dimensional  behavior  than  for  one -dimensional  behavior  to  obtain  the  same  degree  of 
accuracy.  This  is  important  because  practical  complex  structures  such  as  stiffened  plates 
or  shells  are  two-dimensional  and  usually  approximated  by  various  combinations  of  one- 
and  two-dimensional  elements.  Since  the  elements  have  varying  degrees  of  accuracy,  results 
obtained  for  a  structure  approximated  by  a  combination  of  the  elements  may  be  biased  in 
some  sense  rather  than  having  uniform  inaccuracies. 
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SECTION  IV 

CONCLUDING  REMARKS 

Basic  data  are  presented  on  the  convergence  and  accuracy  of  finite  element  equations 
resulting  from  patterns  and  elements  in  common  use.  The  elements  studied  include  bar 
elements,  beam  elements,  plane  stress  elements  of  rectangular  and  triangular  shape,  plate 
bending  elements  of  rectangular  shape  and  straight  and  curved  arch  elements.  The  results 
indicate  that  many  of  the  elements  and  patterns  have  good  convergence  properties;  that  is,  the 
resulting  finite  element  equations  at  a  node  converge  to  the  continuum  equations  at  the  node 
as  the  element  size  vanishes.  The  results  also  indicate  that  some  elements  and/or  patterns 
have  poor  convergence  properties.  Right  triangular  plane  stress  patterns  for  example,  can 
have  undesirable  convergence  properties.  It  is  shown  that  the  requirement  for  interelement 
compatibility  is  neither  necessary  nor  sufficient  to  guarantee  convergence  of  the  finite  element 
equations.  Results  for  arch  approximations  show  that  both  straight  and  curved  elements 
provide  a  convergent  approximation  to  an  arch  structure. 
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SECTION  V 
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SYMBOLS 


A  cross-sectional  area  of  one  dimensional  element 

a  j  b  length  of  rectangular  panel  in  x  and  y  directions,  respectively 


B 

D 

E 

F 


FX,  Fy 


h 

I 

K 

L 

m,  n 


m 


hm  hm 
0(h) 

P 


P  ,  P 
x  y 


q 


R 


t 


u, v,  w 


VVwo 


extensional  stiffness  of  plate,  - 2 

I -ff 

bending  stiffness  of  plate  - - 

)2(l  - /x2) 

Young’s  modulus 
finite  element  nodal  force 

nodal  force  in  x  and  y  directions  of  plane  stress  element 

reference  length  of  finite  element 

moment  of  inertia  of  cross  section 

element  stiffness  matrix 

length  of  one  dimensional  structure 

harmonic  wave  numbers  for  sinusoidally  distributed  loads 

mass  per  unit  length  for  one  dimensional  element;  mass  per  unit  area 
for  two  dimensional  element 

number  of  elements  per  harmonic  wave  number 

omitted  discretization  error  term  of  order  h 

tangential  loading  on  bar  and  arch  structures 

amplitude  of  sinusoidal  in-plane  loading 

distributed  in-plane  loading  in  x  and  y  directions 

transverse  loading 

amplitude  of  sinusoidally  distributed  transverse  loading 
radius  of  arch 
thickness  of  plate 

displacements  in  x,  y  and  z  directions 
amplitude  of  sinusoidally  varying  displacements 
rectangular  coordinates 
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SYMBOLS  (CONT) 


{*} 

€ 

€ 

6,  4> 

H- 

P 


ui 


constant  indicating  ratio  of  element  dimensions 
vector  of  nodal  displacements 

discretization  error  in  displacement  or  frequency  parameter 
discretization  error  in  u  and  v  displacements 
nodal  rotation  variables  for  bending  element  about  y  and  x  axes 
Poisson’s  ratio 

radius  of  gyration  of  cross  section 
circular  frequency 


SUBSCRIPTS 

ex  exact 

it. 

i  identification  for  i  grid  point 
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